home *** CD-ROM | disk | FTP | other *** search
/ c't freeware shareware 1997 / CT_SW_97.ISO / pc / software / wissen / macos / fft121.hqx / FFTs for RISC 1.21 / fftTest.c < prev    next >
C/C++ Source or Header  |  1996-06-27  |  3KB  |  143 lines
  1. /*  A program to test and time complex forward and inverse fast fourier transform routines    */
  2.  
  3. #include <stdio.h>
  4. #include <stdlib.h>
  5. #include <math.h>
  6. #include "fftlib.h"
  7.  
  8. #define    MACTIMING 1            /* true means time fft with macintosh calls    */
  9. #if MACTIMING
  10. #include <timer.h>
  11. #endif
  12.  
  13. #define    NUMROWS 1            /* process Matrix of NumRows different ffts of length N    */
  14. #define N 1024                /* size of FFT must be a power of 2 */
  15. #define NTIMES 20            /* number of timings,invalid if too big (if a[0][0].Re == 0 | nan)*/
  16.  
  17. typedef  struct{
  18.     float Re;
  19.     float Im;
  20.     } Complex;
  21.  
  22. void main(){
  23. float        *Utbl;
  24. Complex        (*a)[N];
  25. #if MACTIMING
  26. UnsignedWide         TheTime1;
  27. UnsignedWide         TheTime2;
  28. UnsignedWide         TheTime3;
  29. double        TheTime;
  30. #endif
  31.  
  32. long         i, il;
  33. long         TheErr;
  34. long        M;
  35.  
  36. Utbl = (float *) calloc((N/4+1),sizeof(float));
  37. if (Utbl==0)
  38.     TheErr = 2;
  39. else
  40. TheErr = FFTInit(&M, N, Utbl);
  41.  
  42. if(!TheErr){
  43.     a = (Complex (*)[N]) calloc(NUMROWS*N,sizeof(Complex));
  44.     if (a == 0) TheErr = 2;
  45. }
  46.  
  47. if(!TheErr){
  48.  
  49.             /*  set up a simple test case */
  50.     for (il=0; il<NUMROWS; il++){
  51.         for (i=0; i<N; i++){
  52.             a[il][i].Re = sqrt(il+i+.77777);    
  53.             a[il][i].Im = (il+i+.22222)*(il+i+.22222) / N - N/2;    
  54.         }
  55.         a[il][0].Re = N+3;
  56.         a[il][1].Re = 1-N;
  57.     }
  58.  
  59. #if MACTIMING
  60.     Microseconds(&TheTime1);
  61.  
  62.     for (i=0;i<NTIMES;i++){        /* do NTIMES times for timing */
  63.         ffts((float *)a, M, NUMROWS, Utbl);
  64.     }    
  65.  
  66.     Microseconds(&TheTime2);
  67.  
  68.     for (i=0;i<NTIMES;i++){        /* do NTIMES times for timing */
  69.         iffts((float *)a, M, NUMROWS, Utbl);
  70.     }    
  71.  
  72.     Microseconds(&TheTime3);
  73.  
  74.  
  75.     TheTime = (double)(TheTime2.hi - TheTime1.hi) * 65536.0 * 65536.0;
  76.     TheTime = (TheTime + (double)(TheTime2.lo - TheTime1.lo))/1000.0;
  77.     printf("Time fft = %12f  ms.", TheTime/NTIMES/NUMROWS, a[0][0].Re);
  78.     TheTime = (double)(TheTime3.hi - TheTime2.hi) * 65536.0 * 65536.0;
  79.     TheTime = (TheTime + (double)(TheTime3.lo - TheTime2.lo))/1000.0;
  80.     printf("  ifft = %12f  ms. a[0][0].Re= %6e\n", TheTime/NTIMES/NUMROWS, a[0][0].Re);
  81.  
  82. #else
  83.     printf("start timing \n");
  84.     for (i=0;i<NTIMES;i++){        /* do NTIMES times for timing */
  85.         ffts((float *)a, M, NUMROWS, Utbl);
  86.         iffts((float *)a, M, NUMROWS, Utbl);
  87.     }    
  88.     printf("end timing \n");
  89. #endif
  90.  
  91.     printf("\n");
  92.             /*  set up a simple test case */
  93.     for (il=0; il<NUMROWS; il++){
  94.         for (i=0; i<N; i++){
  95.             a[il][i].Re = sqrt(il+i+.77777);    
  96.             a[il][i].Im = (il+i+.22222)*(il+i+.22222) / N - N/2;    
  97.         }
  98.         a[il][0].Re = N+3;
  99.         a[il][1].Re = 1-N;
  100.     }
  101.  
  102.     ffts((float *)a, M, NUMROWS, Utbl);
  103.  
  104.     if (N*NUMROWS <= 256){
  105.          for (il=0; il<NUMROWS; il++){
  106.             printf("atrans = [ \n");
  107.             for (i=0; i<N; i++)
  108.                     printf(" %+20.15e + j * ( %+20.15e ) \n", a[il][i].Re, a[il][i].Im);
  109.             printf("]; \n");    
  110.             }
  111.         }
  112.     else { /* abbreviate big output */
  113.         printf("the first fft's last 32 values are: \n");
  114.         for (i=N-32; i<N; i++)
  115.                     printf(" %+20.15e + j * ( %+20.15e ) \n", a[0][i].Re, a[0][i].Im);
  116.     }
  117.  
  118.     iffts((float *)a, M, NUMROWS, Utbl);
  119.  
  120.     if (N*NUMROWS <= 256){
  121.          for (il=0; il<NUMROWS; il++){
  122.             printf("\n aitrans = [ \n");
  123.             for (i=0; i<N; i++)
  124.                     printf(" %+20.15e + j * ( %+20.15e ) \n", a[il][i].Re, a[il][i].Im);
  125.             printf("]; \n");    
  126.             }
  127.         }
  128.     else { /* abbreviate big output */
  129.         printf("\n the first ifft's last 32 values are: \n");
  130.         for (i=N-32; i<N; i++)
  131.                     printf(" %+20.15e + j * ( %+20.15e ) \n", a[0][i].Re, a[0][i].Im);
  132.     }
  133.  
  134.     free (a);
  135.     free (Utbl);
  136.     return;
  137. }
  138. else
  139. if(TheErr==2)    printf(" out of memory ");
  140. else    printf(" error ");
  141. return;
  142. }
  143.